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The differential equations developed by Galejs and How in the analysis of ELF wave 
propagation for frequencies between 60 to 3000 c/s are modified to permit computations 
below the frequency of ion gyrofrequency. Near 8 c/s west-easl propagation exhibits a 
in in in i ii mi ;ilt en nation of 0.022 db/1000 km, while the east-west direction has a local maximum 
of 0.165 db/1000 km. The calculations in the frequency range between 20 and3000c/s are 
made with neglected ion effects, and the propagation is characterized by a monotonic 
decrease of phase velocity down to the lowest frequency. 

The method used for calculating the surface impedance with an arbitrary profile of the 
ionospheric tensor conduct ivity and for strictly east-west , west-east, or south-north direct ions 
of propagation applies also to the VLF range. The VLF propagation characteristics are 
determined using t he modal equat ion for a spherical earth-to-ionosphere geometry. Compu- 
tations based on a simple exponential approximation to the ionospheric conductivity profile 
differ only slightly from computations made using more exact ionosphere parameters. 



1. Introduction 

The propagation of ELF waves below an inhomogeneous anisotropic ionosphere has been 
considered by Galejs and Row [1964|. The propagation is either transverse (east-west and 
west-east propagation) or longitudinal (south-north and north-south propagation) with respeel 
to the horizontal magnetic field. For an arbitrary height variation of the tensor conduc- 
tivity the surface impedance below the horizontally si ratified noncurved ionosphere is com- 
puted by numerically integrating a differential equation. This work has considered the 
frequency range of 60 to 3000 c/s, where the lower frequency limit was determined by the 
accuracy of the numerical integration process. 

In this note the differential equations of Galejs and Row [1964] will be modified in order 
to extend the validity of their numerical solutions to frequencies below the ion gyrofrequency 
(or into the frequency range of earth to ionosphere cavity resonance). The solutions of the 
differential equations also remain valid for higher frequencies and propagation characteristics 
of VLF waves can be determined using an appropriate modal equation. This would extend 
the work of Wait [1963c and 1963d] for isotropic and anisotropic ionosphere of exponential 
conductivity profiles to more general ionosphere models. Other work which has considered 
homogeneous ionosphere models in the VLF range [Wait, 1962, 1963a and b] may provide 
check cases for the present computations. 

Section 2 of this note reviews the relevant ionosphere parameters and sections 3 and 4 
summarize the surface impedance calculations, which follow closely Galejs and Row [1964]. 
The procedure used in calculating the propagation characteristics is out lined in section 5, and the 
numerical results are discussed in section 6. 

2. The Ionosphere Model 

Available data on electron and ion densities and collision frequencies have been sum- 
marized in figures 1 and 2 of Galejs and Row [1964]. These data can be used in the Appleton- 
Hartree formulas for the electrical conductivity of a lightly ionized species. Assuming a 
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Figure 1. Components of the daytime tensor 
conductivity F = 0. 
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Figure 2. Components of the nighttime tensor 
conductivity F = 0. 



static magnetic field B to act in the z direction the Cartesian form of the tensor conductivity 
for a harmonic field with time dependence exp (—iu>t) is, 
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is the plasma frequency of the kih charges species and n k) q k) m k are the number density, charge 
and mass of the kih specie particle. 
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Figure 3. Components of the daytime tensor 
conductivity F = 20 c/s. 
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Figure 4. Components of the nighttime tensor 
conductivity F = 2() c/s. 



is the Arth species gyrofrequency, v k is the effective collision frequency between kth species 
particle and all other constituents of the plasma. Considering electrons (fc=l) in a static 
field of 0.5 oersted (5o=0.5X10" 4 MKS units), co* = 3.185 ^XIO 3 , and ft=-1.04X10 7 . 

The components of the conductivity tensor have been plotted in figures 1 and 2 for /=0, 
in figures 3 and 4 for/=20 and in figures 5 and 6 for /=15 kc/s for both day and night time 
conditions. For /=0 the conductivity components are real, but the components become 
complex for f^0. For frequencies below the ion gyrofrequency the ionic plus electronic con- 
ductivity of /=0 is a good approximation to the real part of the complex conductivity (f=20 
c/s). For frequencies above the ion gyrofrequency the electronic conductivity (the k=\ 
term of (2) to (4)) approximates the real part of the complex conductivity, particularly so for 
altitudes of less than 100 km. In the numerical calculations of propagation constants which 
make use of the electronic conductivity components for/=0, or of the electronic plus ionic 
conductivity for/=15 kc/s, the inequality (x 2 >>o-i is satisfied for altitudes exceeding 100 km. 
This inequality does not apply to electronic plus ionic conductivity for /=() or 20 c/s. 
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Figure 5. Components of the daytime tensor 
conductivity F = 15 kc/s. 
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Figure 6. Components of the nighttime tensor 
conductivity F=15 kc/s. 



3. Equatorial Propagation 

For an assumed exp (—iwt) time dependence of the fields the relative tensor permittivity 
is of the form 
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if the static magnetic field is in the z direction. It is assumed that the field components 
exhibit no variation in the z direction and that [e] is changing only in the y direction, which 
would be the south-north and vertical directions in terrestrial propagation problems. For 
east to west propagation along the magnetic equator the surface impedance below the aniso- 
tropic medium is given by Galejs and Row [1964] as 
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and 



The sign of e 2 is reversed for west to east propagation, which causes an obvious nonreciprocity. 
For propagation in the x direction the magnetic field component H z will exhibit the x and y 
dependence of 

tf 2 = exp (ik x x)Y(y). (9) 

Y{y) is determined from a solution of the differential equation 
where 

P( ^| lo K^t) (n) 

QM-Z&K-M*^^) (12) 

with k = «V ^oMo- The solution to this equation is formally given by 

) • ( y) = e -lS^f (y)= (__^- U2 f(y) (13) 

where f(y) satisfies the equation 

dy 

with 

™=tf- 2 f-i p2 - (15) 

The surface impedance (S) can be expressed in terms of /'(//) as 

Z s = ^-U^h) I V{y)-\ P(y)~}-^- 2 ik x \ (16) 



/ - / 2 +A'(//) = (14) 



where 



V(y)= l y d {. (17) 



F(t/) satisfies the Riccati-type differential equation 

~-[V(yy+K(y)}. (18) 

However, V(y) becomes approximately equal to \ P(y) for the lower frequencies near the lower 
ionosphere boundary, and minor errors in the numerical integration cause significant inac- 
curacies in Z s . This difficulty can be avoided by defining a new variable 

U{y)=V{y)-\p(y). (19) 

Substitution of (19) in (IS) shows that U(y) satisfies the differential equation 

( ^ / = -U- 2 (y)-P(y)U(y)-Q(y), (20) 

which can be integrated numerically if U=U(y u ) is specified at a given altitude y=y u . U(y u ) 
may be determined from the surface impedance Z s at the upper edge y=y u of the ionosphere 
model. Examination of the conductivity profiles of figures 1 and 2 indicates that o- 2 >>o-i at 
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y>100 km with ions neglected or for o>>6r ion with ions considered. Assuming an exponential 
variation of a 2 for y^>y U ^>100 km and setting e x = constant, U(y u ) is computed as 

U(y u )=ik ^. (21) 

It may be also assumed that ci and c 2 = constant for y^>y u - The corresponding value of 
U{y u ) is 

U{y u )=ihyj^, (22) 

which is applicable to any ratio of o^M- The above equations apply to east to west propaga- 
tion. For west to east propagation e 2 in (12), (15), and (16) should be multiplied by minus one. 

4. Propagation Across the Magnetic Equator 

Near the magnetic equator the static earth's magnetic field is in the south to north direc- 
tion, which is the 2-direction in the definition of the permittivity tensor (7). The tensor [e] 
changes only in the vertical y-direction and the field components exhibit no variation in the 
^-direction (west to east direction). The differential equations of the various field components 
are generally coupled. However, in the lower ionosphere the coupling may be neglected 
[Galejs and Row, 1964] and the surface impedance for vertically polarized fields is simply 

z ° = ETi^K^ Hx - (23) 

H x is of the form 

fl*=exp (ik,z)Y(y). (24) 

Y(y) is determined from the solution of the differential equation 

SlHtf +"«-*>r-* (25) 

The parameter q of (25) characterizes the amount of coupling between TE and TM field com- 
ponents (or between E x and H x ) at a given altitude y. There is negligible coupling with q=l 
near the lower ionosphere boundary, but the amount of coupling is increased and the magnitude 
of q is decreased with increasing altitude. The simplest approximation of 

implies negligible coupling between E x and H x for y<Cyi, arid strong coupling for y>y lt This 
artificial boundary at y=yi will introduce an impedance discontinuity at a height where it 
appears to have little effect on fields below. In analogy with section 3, the surface impedance 
is given by 

Z s =^~\v{y)-\l\y)]=^- l -U(y) (27) 

?,coeoe 3 L 2 J ^coeo€3 

where 

P(2/)=-^(lo ge8 ). (28) 

U(y) is determined from the solution (20) where Q(y) is given by 

Q<3f)=+R--q%. (29) 
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The differential equation (20) is integrated downward from y=y u ^>y 1 . The starting value 
U(y u ) is computed for an exponential height variation of e 3 and for e x = constant. It follows 
that 

U(yJ=H%-qkilei<Jei. (30) 

In view of (26) U(y u ) of (30) is changed to 

U{y u )=ih<JT z . (31) 

5. Propagation Characteristics 

5.1. ELF 

For ELF waves the propagation constant ik x of section 3 or ik z of section 4 is given by 
ik S, where S is computed following Galejs [1961] from an iteration solution of 



where 



s » = V 1+ j^k^)V 1- [^fe] (32) 



Z s is computed at the assumed boundary height y of the ionospheric profile. Z s of (8), (16)? 
(23), or (27) depend also on S. The subscripts of S designate the successive iterations, and the 
calculations may be started with S Q =l in the right-hand side of (32). 

5.2. VLF 

The propagation characteristics of the VLF waves are calculated starting out with 
(VII-27) of Wait [1962], which is rewritten for exp (-— i<d) time dependence as 

r w^-g.wAtn r w2(t-u)+ qi w 2 (t-u) i_^ 2irn 

L^2(t) — q g w 2 (t)j\jv[(t—u) + q i Wi(t--u)J K } 

-t=(*fj K C>, (35) 
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Wi(£) and w 2 (t) are the Airy functions, n is the mode number, a is the radius of the earth, and 
C is a root of the modal equation (6 r2 =l— S 2 ). The Airy functions may be represented by a 
power series or they can be related to Bessel or Hankel functions [Spies and Wait, 1961 1. The 
Airy functions of arguments (t — u) are expressed in terms of Hankel functions as in (VI 1—42) 
of Wait [1962], which are .replaced by asymptotic expansions. This gives 

log (-BO+h* ffglgggj +i § MO*-* \ (4»-l)-0 (40) 
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Approximate solutions of (40) can be worked out on the assumptions k^C'Y^^ 1, (C'/At) 
<0, ^ T2< C< (2y/ a ) and ^<<1. Approximating the term of (40) which involves Airy functions 
by its (7=0 value, expanding log (—Ri) in powers of W'/Ai) and introducing a further 
approximation 
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where M must be such that the approximation (43) is still valid. For (C'/Ai) << 1 only the 
m = terms of (44) are significant and (44) can be seen to become equivalent to (IX-19) of 
Wait [1962]. 

More accurate solutions of (40) are obtained by Newton's method. A Taylor series 
expansion of a function j\z) about its approximate zero J(z ) gives a corrected z estimate for 
its zero as 

/(go) 



z = z — 



dj/dz\ 



(45) 



Iterations may be used to improve accuracy of this process. Letting z=0 2 J the initial value z 
is obtained from (44) by assuming a perfectly reflecting ionosphere of A^oo (or R i—> — 1). The 
function f(z ) and its derivative follow from (40) as 
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For a=0=O and <7«A 4 (46) to (49) are equivalent to (VII-37) to (VII-40) of Wait [1962]. 
The propagation constant ik x of section 3 or ik z of section 4 is given by ik S where 



s=Vi-c 2 . 

700 



(50) 



The surface impedances Z s of (16) and (27) are computed for the assumed lower boundary y of 
the ionospheric profile. Z s is also dependent on S. S of (50) is computed using iterations by 
setting S=\ in the initial Z & expression and by updating the S value in Z, following each £ 
computation. For modes of low attenuation the process is rapidly convergent and only a 
few iterations are necessary. 

The above computational procedure differs from the method used by Wait [1963c] where 
the modal equation, similar to (44), is considered at a height (y-\-5) which is above the lower 
boundary y of the ionosphere. A surface impedance is computed in a planar geometry at y, 
and the reflection coefficient Ri(y) is referred to the height (y-\-8) as for free space in this 
planar increment 8 [Wait and Walters, 1963]. The reflection coefficient Rt(y-\-8) becomes 

R i (y+5)=R l (y) exp (-i2k Cd). (51) 

This is equivalent to transforming Ai(y) = A t into A z (t/ + <5) by 

Substituting (52), M=0 and assuming <5<^<C?/, (44) may be rearranged into 

which will differ from C 2 for 5^0, particularly if 5 is a sizable fraction of ;?/. The method which 
computes the surface impedance at ?/, but considers the modal equation at ?/ + <5 [Wait, 1963c], 
will therefore differ from computations which consider the modal equation also at 7/. 1 The 
latter approach is strictly correct only if the surface 4 impedance is computed in a spherical 
geometry, or if the planar surface impedance is a good approximation to the spherical one. 
Impedance transformations similar to (52) but suitable for a spherical geometry are discussed 
in the appendix. The impedance error due to planar calculations can be decreased by select- 
ing 5 to be different from zero, although no such calculations are reported in this paper. 

6. Discussion 

The roots of the modal equation and the ELF attenuation rates have been calculated for 
the ionospheric conductivity profiles of figure 1 to 4. The ionospheric profile with considered 
ions for f=20 c/s applies to propagation below the ion gyrofrequency and was used for calculat- 
ing the curves of figures 7 and 8. The plot in figure 7 shows a minimum of daytime Im S for 
east-west propagation around 8 c/s. Im S for west-east direction has a maximum at this 
frequency. The corresponding attenuation rate figures are a^ = 0.165 db/1000 km, and a we = 
0.022 db/1000 km. The field enhancement near the first resonant frequency of the earth-to- 
ionosphere cavity appears, therefore, due more to west-east than to east-west directions of 
propagation. The roots of the modal equations and the attenuation rates have been also 
calculated for electronic conductivity (f=0) in the frequency range of 20 to 2000 c/s and are 
shown in figures 9 to 11. It was already pointed out that this appears to be a satisfactory 
approximation to the computed conductivity profiles for frequencies above the ion gyro- 
frequency. In distinction from earlier work [Galejs and Row, 1964], the Re S is now mono- 



i Such col 1 1 it mmi is do not apply to planar geometry (e.g., eq VI-4.1 of Wait [1962]) where a change of the boundary height of the modal equation 
and a simultaneous transformation of Ri as in (51) do not affect the solution. However, the planar mode equation is not considered accurate in 
the VLF range. 
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Figure 7. Roots of the modal equation. 
Daytime conductivity with ions considered (fig. 3). 
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Figure 9. Roots of the modal equation. 
Electronic daytime conductivity for f=0 (fig. 1). 



1 1 1 


' 1 1 1 I 




1 




- 


t 


mmm — — — Im 5 


- 


\ 




- 


\ 






\ 






\ 






\ 






. EW 






\ 




_ 


VWE \ 




- 




"V 




^NS ^^h^ 




- 




r-^^-^^^j^-- ._ ___ 




WE 




*^- ^ — i 




NS 


*^ "-^^!* 


_ 


/^EW 




- 


1 1 1 


i 1 i 1 1 





20 30 

FREQUENCY - C/S 



Figure 8. Roots of the modal equation. 
Nighttime conductivity with ions considered (fig. 4). 
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Figure 10. Roots of the modal equation. 
Electronic nighttime conductivity for f=0 (fig. 2). 



tonically increasing or the phase velocity monotonically decreasing down to the lowest fre- 
quencies of the plot. The high-frequency behavior of these curves is about the same as in the 
computations by Galejs and Row [1964]. 

The calculations in the VLF range were first made for a sharply bounded homogeneous 
ionosphere to obtain a check with published results. The ionospheric conductivity is repre- 
sented by the tensor (7). The computed curves are similar to those of Wait [1962], which 
were obtained by a different method. East-west direction exhibits the higher attenuation. 
The differences between the Re S figures are almost indistinguishable, although Re S ew exhibits 
the higher values. The curves for an isotropic ionosphere are intermediate between those of 
east-west and west-east directions. 
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Figure 12. Re $ for exponential ionosphere model. 
It is assumed that free space exists below the indicated values of y. 

Similar remarks apply to figures 13 and 14. 
Figure 11. Attenuation constants. Electronic con- 

ductivity for f=0 (figs. 1 and 2). 

The height of the lower ionosphere boundary which was used for surface impedance calcu- 
lations is designated by ?/-consL in the figures. 

The exponential ionosphere is assumed to have a conductivity profile of 

(r,=10- 6 exp |().291(?/ -74)] (54) 

where j=l or 2 in the definition of the permittivity tensor (7) and y is in kilometers. The 
effects of frequency changes are shown in figures 12 and 13. Increasing the height of the lower 
ionosphere boundaries increases Re S or decreases the phase velocity (v p /c=l/Re S) similarly 
as for the sharply bounded ionosphere. The attenuation rates show a significant change with 
different assumed heights of the lower ionosphere boundaries for the same ionospheric profile 
(54). These attenuation rates, obtained from the solution of (45) to (48) are compared in 
figure 14 with the approximate solution (44). There is a qualitative agreement, although 
differences of approximately 0.5 d\) may be noted. 

The effects of profile shifts are illustrated in figures 15 and 16. An upward shift of the 
profile (<5>0) increases Re S and decreases attenuation for a constant boundary height. The 
same conclusion applies in most cases to simultaneous profile and boundary shifts (e.g., 5 = 
and +5 curves are compared for boundaries at y and y+5 respectively). However, increasing 
the boundary height of a profile increases the attenuation rates significantly for y approaching 
80 km. 

These characteristics of attenuation rates apply to east-west, west-east propagation as 
well as to the isotropic ionosphere, and will be discussed with the aid of a surface impedance 
curve for the isotropic model shown in figure 17. The numerically calculated surface imped- 
ance can be checked against a closed form solution [Wait 1962, Galejs 1961] which is strictly 
valid for large magnitude of the local refractive index (\n\ = |l+io-/we | 1/2 >>l). The closed 
form solution agrees with numerical calculations for ?y> 80 km, where the earlier nonequality 
is satisfied. For ?/<60 km the computer calculation and the free space approximation (52) 
both give about the same y dependence of A*. The normalized surface impedance |A ?: | is 
decreased in the height interval of 75 to 80 km where Re A t « |Im A*|. A decreasing Re A t 
gives an increase of (— Im C 2 ) and hence of (Im 8) or a following (44) which accounts for 
the a increase shown in figure 16. There is a minimum of a near y=70 km where Re A< is 
maximum in figure 17. This reasoning does not apply to ?/<60 where Re A/<lm A*, 
and where small Re A* does not imply a large Im S or a. 
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''igure 13. Attenuation rates for exponential iono- 
sphere model. 
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Figure 15. Re S for shifted exponential profiles. 
Free space conditions are assumed below the boundary height. 
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Figure 14. Attenuation rates for exponential iono- 
sphere model: comparison of approximations with 
more exact solutions. 
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Figure 16. Attenuation rates for shifted exponential 

profiles. 
Free space conditions are assumed below the boundary height. 



Further calculations have been made for the daytime ionosphere conductivity profile of 
figure 5. The general behavior of the propagation parameters in figures 18 and 19 is similar to 
that shown in figures 12 and 1 3 for the simple exponential profile. A comparison with the approxi- 
mate solution (44) is shown in figure 20, and comparisons with measurements [Watt and 
Croghan, 19631 are indicated in figure 21. There is a qualitative agreement although calcula- 
tions give a nearly 1 db lower attenuation rate for north-south and west-south directions of 
propagation. 
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Figure 17. Surface impedance for the isotropic 
exponential ionosphere model. 
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Figure 19. Daytime attenuation rates; electronic plus 
ionic conductivity (fig. 5). 
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Figure 18. Re S for daytime ionosphere; electronic 

plus ionic conductivity (fig. 5) . 

It is assumed that free space exists below the indicated values of y. 

Similar remarks apply to figures 10, 20, and 21 . 
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y =60 = CONST. 
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Figure 20. Daytime attenuaticn rates; electronic 
plus ionic conductivity (fig. 5). 
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THEORY (y = 60 km) 



MEASUREMENTS 




FREQUENCY-KILOCYCLES 



Figure 21. Measured and calculated daytime 
attenuation rates. 



Calculations based on nighttime ionosphere profiles of figures 2 or 6 are less successful. 
Isotropic propagation results in a approximately 2 db/1000 km, the west-east and east-west 
directions exhibit a<0.5 and a>3 db/1000 km respectively. It is particularly difficult to 
justify the low a figure for the west-east direction of propagation and these numerical results 
are not shown in detail. 

This work on the VLF propagation represents an initial attempt to consider a numerical 
surface impedance calculation in conjunction with the modal equations for spherical earth. 
The simple exponential conductivity profile gives about the same propagation constants as 
more elaborate ionosphere models, but the calculations are in both cases critically dependent 
on the assumed lower boundary height of the ionosphere model. A possible cause for this is 
the probable incompatibility between planar surface impedance calculations and the particular 
modal equation for spherical geometry. However, some improvement may be expected by 
referring the calculated planar surface impedance to a height near the equivalent reflection 
level of the ionosphere [Wait, 1963c]. It appears that further work is needed to clarify fully 
the effects of various boundary heights in the models of nonhomogeneous anisotropic ionosphere. 



Appreciation is expressed to E. M. Larsen and J. Osborn for computational assistance. 
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7. Appendix 

Surface Impedance Transformation in a Spherical Geometry 

The surface impedance computed at the bottom of the ionosphere y may be referred to a 
higher altitude (y+8) if this does not alter the solutions of the modal equation. Such surface 
impedance transformations will be discussed in this appendix for the modal equations of 
section 4.2. 
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The modal equation (40) with /?, as in (41), but with neglected terms proportional to 
(C) -3 has constant solutions if 



At-C 



', exp [ ^ (0") 3 ]=const, (55) 



At+C 
This is equivalent to transforming A t (y) = Ai to A, ; (?/ + <5) by 

i i+crffi^Mcm , (y) 



A t (y+8) 1+ iCM tm[koSO , iy)] C>(y+S) 



(56) 



with 

C'(y) 8( C*a\ 

C'iy+sr 1 2y{ 1 y ) (57) 

For £„5( "(?/)< <1 (56) simplifies to 

z^=a^)[-f?('+f)-|('-f)] 

However, this form of impedance transformation does not provide a constant solution of the 
more approximate modal equation (44). The modal equation (44) with M=0 remains un- 
altered if 

^—\ iha(Cy=con$t. (59) 



This gives 



i^=(}+^)<^ v» 



where C'(y)/C'(y+8) is given by (57). For a^oo, r;'(?/)->C r/ (?/-f<5) following (42), and (60) 
reduces to (52). When applying- (60) and (57) to (44), it is essential to keep the term propor- 
tional to 6 y2 in (57). 

These surface impedance transformations should be used in connection with surface 
impedance calculations for a spherically stratified medium. 
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